function [DiagVec, DiagBlock] = ConstructHMatrix(Kernel, nLevel)

global RKRow RKCol RKBlock nList;

nx = length(Kernel);
nWidthTop = nx / 4;
% nWidthBottom = 32;
% nLevel = round( log2( nWidthTop / nWidthBottom ) + 1);

maxRank = 4;
nSample = maxRank + 20;


for iLevel = 1 : nLevel
  % Setup. The top level is different
  if( iLevel == 1 )
    nWidthCur = nWidthTop;
    nBlockCur = 4;
    nRKBlockPrevSum = 0;
    nList = 2;
    nRKBlock = 2;
    RKRow = zeros( nRKBlock, 2 );
    RKCol = zeros( nRKBlock, 2 );
    RKBlock = cell( nRKBlock, 1 );
   
    RKRow(1, :) = [1 nWidthCur];
    RKRow(2, :) = [nWidthCur+1 2*nWidthCur];
    RKCol(1, :) = [2*nWidthCur+1 3*nWidthCur];
    RKCol(2, :) = [3*nWidthCur+1 4*nWidthCur];
  
    for iRK = 1 : nRKBlock
      RKBlock{iRK} = struct('U', [], 'V', []);
      RKBlock{iRK}.U = zeros( nWidthCur, maxRank );
      RKBlock{iRK}.V = zeros( nWidthCur, maxRank );
    end

    SampleList = cell( nList, 1 );
    for iList = 1 : nList 
      SampleList{iList} = iList;
    end
  else
    nWidthCur = nWidthCur / 2;
    nBlockCur = nBlockCur * 2;
    nRKBlockPrevSum = nRKBlockPrevSum + nRKBlock
    nList = 12;
    nRKBlock = nList * 2^(iLevel-2);
    
    RKRow = [RKRow; zeros( nRKBlock, 2 )];
    RKCol = [RKCol; zeros( nRKBlock, 2 )];
    RKBlock = [RKBlock; cell( nRKBlock, 1 )];

    SampleList = cell( nList, 1 );

    % This could be improved later
    for iList = 1 : nList 
      SampleList{iList} = nRKBlockPrevSum + iList : nList : ...
	nRKBlockPrevSum + nRKBlock;
    end


    c = nRKBlockPrevSum;
    for iRow = 1 : nBlockCur
      if( mod( iRow, 2 ) == 1 ) 
	c = c + 1;
	RKRow(c, :) = [(iRow-1)*nWidthCur + 1 iRow*nWidthCur]; 
	c = c + 1;
	RKRow(c, :) = RKRow(c-1, :);
      else
	c = c + 1;
	RKRow(c, :) = [(iRow-1)*nWidthCur + 1 iRow*nWidthCur];
      end
    end

    c = nRKBlockPrevSum;
    for iRow = 1 : nBlockCur
      ind = mod( iRow + 1, nBlockCur ) + 1;

      if( mod( iRow, 2 ) == 0 ) 
	c = c + 1;
	RKCol(c, :) = [(ind-1)*nWidthCur+1  (ind)*nWidthCur];
	c = c + 1;
	RKCol(c, :) = RKCol(c-1, :);
      else
	c = c + 1;
	RKCol(c, :) = [(ind-1)*nWidthCur+1  (ind)*nWidthCur];
      end
    end

  end % if; setup

  % Computation part
  for iList = 1 : nList
    TempR = zeros( nWidthCur, nSample, length(SampleList{iList}));
    TempQ = zeros( nWidthCur, maxRank, length(SampleList{iList}));
    TempT = zeros( nWidthCur, maxRank, length(SampleList{iList}));

    for iSample = 1 : nSample
      
      
      % Sample
      vecRand = zeros( nx, 1 );
      for iRK= 1 : length(SampleList{iList})
        vecRand( RKRow( SampleList{iList}(iRK), 1 ) : RKRow( ...
            SampleList{iList}(iRK), 2 )) = randn( nWidthCur, 1 );
      end
      vecResult = KernelApply( Kernel, vecRand );

      % Subtract the low rank matrices from previous levels. This step is not
      % necessary for the second level.

      vecResultRK = SuperMatrixApply( vecRand, nRKBlockPrevSum);
      vecResult = vecResult - vecResultRK;	

      % Distribute vecResult into interaction blocks
      for iRK = 1 : length(SampleList{iList})
        TempR(:, iSample, iRK) = vecResult( ...
            RKCol( SampleList{iList}(iRK), 1 ) : ...
            RKCol(SampleList{iList}(iRK), 2 ) );
      end % iRK
    end %iSample

    for iRK = 1 : length(SampleList{iList})
      [TempQ1, TempS1, TempV1] = svd(conj(TempR(:, :, iRK)), 'econ');
      TempQ(:, :, iRK) = TempQ1(:, 1:maxRank );
    end % iRK
    clear TempR;

    for iRank = 1 : maxRank
      vecQ = zeros(nx, 1);
      for iRK = 1 : length(SampleList{iList})
        vecQ( RKCol( SampleList{iList}(iRK), 1 ) : RKCol( ...
            SampleList{iList}(iRK), 2 ) ) = TempQ(:, iRank, iRK);
      end
      vecResult = KernelApply( Kernel, vecQ );

      vecResultRK = SuperMatrixApply( vecQ, nRKBlockPrevSum );
      vecResult = vecResult - vecResultRK;


      % Distribute vecResult into interaction blocks

      for iRK = 1 : length(SampleList{iList})
        TempT(:, iRank, iRK) = vecResult( ...
            RKRow( SampleList{iList}(iRK), 1 ) : RKRow( ...
                SampleList{iList}(iRK), 2 ) );
      end %iCol
    end	% iRank

    for iRK = 1 : length(SampleList{iList})
      [TempU1, TempS1, TempW1 ] = svd( TempT(:, :, iRK), 'econ' );
      TempV1 = TempQ(:, :, iRK) * TempW1 * TempS1;
      
      RKBlock{SampleList{iList}(iRK)}.U = TempU1;
      RKBlock{SampleList{iList}(iRK)}.V = TempV1;
      
    end %iRK
    
    if (0)
      RKRow(SampleList{iList}(iRK), :)
      RKCol(SampleList{iList}(iRK), :)
      KernelExact = Kernel( ...
          RKRow(SampleList{iList}(iRK), 1) : RKRow(SampleList{iList}(iRK), ...
                                                   2), ...
          RKCol(SampleList{iList}(iRK), 1) : RKCol(SampleList{iList}(iRK), ...
                                                   2));
      K = TempU1 * (TempV1)';
      norm(K - KernelExact)
%      S  = svd(KernelExact);
%      semilogy(S(:))
    end
    
    clear TempQ, TempT; 
    
%    pause

  end % iList

end


% Extract diag part

nRKBlockPrevSum = nRKBlockPrevSum + nRKBlock;

% nWidthCur = nWidthCur;
nList = nWidthCur * 2;
nSample = nx / nList;
DiagBlock = cell( nSample*2, 1 );
for i = 1 : 2*nSample
  DiagBlock{i} = zeros( nWidthCur, nWidthCur );
end
DiagVec = zeros( nx, 1 );

for iList = 1 : nList
  vecCol = zeros(nx, 1);
  vecCol(iList:nList:end) = 1;
  vecResult = KernelApply( Kernel, vecCol );
  vecResultRK = SuperMatrixApply( vecCol, nRKBlockPrevSum );
 
  % figure(1)
  % clf
  % hold on
  % plot(real(vecResult),'b')
  % plot(real(vecResultRK),'r')
  % xlim([200 250])
  % hold off
  % pause
  % vecResultRK = zeros( nx, 1 );
  
  vecResult = vecResult - vecResultRK;
  DiagVec(iList:nList:end) = vecResult(iList:nList:end);
  if( iList <= nWidthCur )
    for i = 1 : nSample
      DiagBlock{2*i-1}(:, iList) = ...
	vecResult( 1+((2*i-1)-1)*nWidthCur: (2*i-1) * nWidthCur );
    end 
  else
    for i = 1 : nSample
      DiagBlock{2*i}(:, iList-nWidthCur) = ...
	vecResult( 1+((2*i)-1)*nWidthCur: (2*i) * nWidthCur );
    end 
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = KernelApply( A, u )
v = A*u;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = SuperMatrixApply( u, nRKBlockPrevSum)
global RKRow RKCol RKBlock;
eps = 1d-15;
nx = length(u);
v = zeros( nx, 1 );
for iList = 1 : nRKBlockPrevSum
  if( norm( u( RKCol(iList, 1): RKCol(iList,2)) ) > eps )
    v(RKRow(iList, 1): RKRow(iList,2)) = ...
      v(RKRow(iList, 1): RKRow(iList,2)) + ...
      RKBlock{iList}.U * ...
      ((RKBlock{iList}.V)' * u( ...
      RKCol(iList, 1): RKCol(iList,2)) );
  end
  % Symmetrize
  if( norm( u( RKRow(iList, 1): RKRow(iList,2)) ) > eps )
    v(RKCol(iList, 1): RKCol(iList,2)) = ...
      v(RKCol(iList, 1): RKCol(iList,2)) + ...
      conj(RKBlock{iList}.V) *  ...
      (transpose(RKBlock{iList}.U) * u( ...
      RKRow(iList, 1): RKRow(iList,2)) );
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

